program brentry_bioprobit_neg_3t
version 10.1
args lnf s1 v1 f1 alpha1_1 gamma1_1 s2 v2 f2 alpha1_2 alpha2_2 alpha3_2 gamma1_2 gamma2_2 gamma3_2 s3 v3 f3 alpha1_3 alpha2_3 gamma1_3 gamma2_3 r1 r2 r3

tempvar pi_h1 pi_a1 pi_a2 pi_u1 pi_u2 pi_u3 p1u p2u p3u p1h p1a p2a p110 p120 p130 p101 p102 p011 p021 p031 p012 p022 p032 p111 p121 p131 p112 p122 p132 
tempname rho1 rho2 rho3 Sigma

qui gen double `rho1' = tanh(`r1')
qui gen double `rho2' = tanh(`r2')
qui gen double `rho3' = tanh(`r3')

qui gen double `pi_h1' = `s1'*(`v1' + `alpha1_1') -`f1' - `gamma1_1'

qui gen double `pi_a1' = `s3'*(`v3' + `alpha1_3') - `f3' - `gamma1_3'
qui gen double `pi_a2' = `s3'*(`v3' + `alpha1_3' - `alpha2_3') - `f3' - `gamma1_3' - `gamma2_3'

qui gen double `pi_u1' = `s2'*(`v2' + `alpha1_2') -`f2' - `gamma1_2'
qui gen double `pi_u2' = `s2'*(`v2' + `alpha1_2'-`alpha2_2') -`f2' - `gamma1_2' - `gamma2_2'
qui gen double `pi_u3' = `s2'*(`v2' + `alpha1_2'-`alpha2_2'-`alpha3_2') -`f2' - `gamma1_2' - `gamma2_2' - `gamma3_2'

qui gen double `p1u'=normal(-1*`pi_u1')
qui gen double `p2u'=normal(-1*`pi_u2')
qui gen double `p3u'=normal(-1*`pi_u3')

qui gen double `p1h'=normal(-1*`pi_h1')
	
qui gen double `p1a'=normal(-1*`pi_a1')	
qui gen double `p2a'=normal(-1*`pi_a2')	

qui gen double `p110'=binormal(-1*`pi_h1', -1*`pi_u1', `rho1')
qui gen double `p120'=binormal(-1*`pi_h1', -1*`pi_u2', `rho1')
qui gen double `p130'=binormal(-1*`pi_h1', -1*`pi_u3', `rho1')

qui gen double `p101'=binormal(-1*`pi_h1', -1*`pi_a1', `rho2')
qui gen double `p102'=binormal(-1*`pi_h1', -1*`pi_a2', `rho2')

qui gen double `p011'=binormal(-1*`pi_u1', -1*`pi_a1', `rho3')
qui gen double `p021'=binormal(-1*`pi_u2', -1*`pi_a1', `rho3')
qui gen double `p031'=binormal(-1*`pi_u3', -1*`pi_a1', `rho3')

qui gen double `p012'=binormal(-1*`pi_u1', -1*`pi_a2', `rho3')
qui gen double `p022'=binormal(-1*`pi_u2', -1*`pi_a2', `rho3')
qui gen double `p032'=binormal(-1*`pi_u3', -1*`pi_a2', `rho3')

local m_rho1 = tanh(`r1')
local m_rho2 = tanh(`r2')
local m_rho3 = tanh(`r3')
qui matrix `Sigma' = (1, `m_rho1', `m_rho2' \ `m_rho1', 1, `m_rho3' \ `m_rho2', `m_rho3', 1)

qui gen double `p111' = .
mata: my_mvnorm("`pi_h1'", "`pi_u1'", "`pi_a1'", "`p111'", "`Sigma'")
qui gen double `p121' = .
mata: my_mvnorm("`pi_h1'", "`pi_u2'", "`pi_a1'", "`p121'", "`Sigma'")
qui gen double `p131' = .
mata: my_mvnorm("`pi_h1'", "`pi_u3'", "`pi_a1'", "`p131'", "`Sigma'")
qui gen double `p112' = .
mata: my_mvnorm("`pi_h1'", "`pi_u1'", "`pi_a2'", "`p112'", "`Sigma'")
qui gen double `p122' = .
mata: my_mvnorm("`pi_h1'", "`pi_u2'", "`pi_a2'", "`p122'", "`Sigma'")
qui gen double `p132' = .
mata: my_mvnorm("`pi_h1'", "`pi_u3'", "`pi_a2'", "`p132'", "`Sigma'")

qui replace `lnf' = ln(`p111') if $ML_y1==0 & $ML_y4==0 & $ML_y7==0
qui replace `lnf' = ln(`p121'-`p111') if $ML_y1==0 & $ML_y4==1 & $ML_y7==0
qui replace `lnf' = ln(`p131'-`p121') if $ML_y1==0 & $ML_y4==2 & $ML_y7==0
qui replace `lnf' = ln(`p101'-`p131') if $ML_y1==0 & $ML_y4>=3 & $ML_y7==0

qui replace `lnf' = ln(`p011'-`p111') if $ML_y1>=1 & $ML_y4==0 & $ML_y7==0
qui replace `lnf' = ln(`p021'-`p011'-`p121'+`p111') if $ML_y1>=1 & $ML_y4==1 & $ML_y7==0
qui replace `lnf' = ln(`p031'-`p021'-`p131'+`p121') if $ML_y1>=1 & $ML_y4==2 & $ML_y7==0
qui replace `lnf' = ln(`p1a'-`p031'-`p101'+`p131') if $ML_y1>=1 & $ML_y4>=3 & $ML_y7==0

qui replace `lnf' = ln(`p112'-`p111') if $ML_y1==0 & $ML_y4==0 & $ML_y7==1
qui replace `lnf' = ln(`p122'-`p121'-`p112'+`p111') if $ML_y1==0 & $ML_y4==1 & $ML_y7==1
qui replace `lnf' = ln(`p132'-`p131'-`p122'+`p121') if $ML_y1==0 & $ML_y4==2 & $ML_y7==1
qui replace `lnf' = ln(`p102'-`p132'-`p101'+`p131') if $ML_y1==0 & $ML_y4>=3 & $ML_y7==1

qui replace `lnf' = ln(`p012'-`p011'-`p112'+`p111') if $ML_y1>=1 & $ML_y4==0 & $ML_y7==1
qui replace `lnf' = ln(`p022'-`p021'-`p012'-`p122'+`p011'+`p121'+`p112'-`p111') if $ML_y1>=1 & $ML_y4==1 & $ML_y7==1
qui replace `lnf' = ln(`p032'-`p031'-`p022'-`p132'+`p021'+`p131'+`p122'-`p121') if $ML_y1>=1 & $ML_y4==2 & $ML_y7==1
qui replace `lnf' = ln(`p2a'-`p1a'-`p032'-`p102'+`p031'+`p101'+`p132'-`p131') if $ML_y1>=1 & $ML_y4>=3 & $ML_y7==1

qui replace `lnf' = ln(`p110'-`p112') if $ML_y1==0 & $ML_y4==0 & $ML_y7>=2
qui replace `lnf' = ln(`p120'-`p110'-`p122'+`p112') if $ML_y1==0 & $ML_y4==1 & $ML_y7>=2
qui replace `lnf' = ln(`p130'-`p120'-`p132'+`p122') if $ML_y1==0 & $ML_y4==2 & $ML_y7>=2
qui replace `lnf' = ln(`p1h'-`p130'-`p102'+`p132') if $ML_y1==0 & $ML_y4>=3 & $ML_y7>=2

qui replace `lnf' = ln(`p1u'-`p110'-`p012'+`p112') if $ML_y1>=1 & $ML_y4==0 & $ML_y7>=2
qui replace `lnf' = ln(`p2u'-`p1u'-`p120'+`p110'-`p022'+`p012'+`p122'-`p112') if $ML_y1>=1 & $ML_y4==1 & $ML_y7>=2
qui replace `lnf' = ln(`p3u'-`p2u'-`p130'+`p120'-`p032'+`p022'+`p132'-`p122') if $ML_y1>=1 & $ML_y4==2 & $ML_y7>=2
qui replace `lnf' = ln(1-`p1h'-`p3u'+`p130'-`p2a'+`p032'+`p102'-`p132') if $ML_y1>=1 & $ML_y4>=3 & $ML_y7>=2

end

mata
void my_mvnorm(string scalar hname, string scalar uname, string scalar aname, string scalar pmvname, string scalar sigma)
{
	Y = vech(st_matrix(sigma))'
	R = J(st_nobs(), 1, Y)
	x = st_data(., hname)
	y = st_data(., uname)
	z = st_data(., aname)
	st_view(pmv=., ., pmvname)
	U = -1*x, -1*y, -1*z
	pmv[.,.] = mvnormal(U, R)
}
end